###########################################
### SUPPLEMENTALS                       ###
###########################################

# Data preparation --------------------------------------------------------

# load packages

pckgs <-
  c(
    "raster",
    "grid",
    "gridExtra",
    "dplyr",
    "ggplot2",
    "devtools",
    "Hmisc",
    "glmmTMB",
    "sjmisc",
    "sjlabelled",
    "sjPlot",
    "sjstats",
    "ggeffects",
    "parameters",
    "stargazer")

lapply(pckgs, require, character.only = TRUE)

# set wd

setwd("C:/Users/majacob/polybox/Statistics/R/2021_party_income_states/replication")

# loading data

data  <-
  read.csv(file = "Germany_finance_states_inc_analysis_2009-2017.csv", na.strings =
             c("", "NA"))

# Figure A1 ---------------------------------------------------------------

# Compute share

data$total_income <- data$i_fees+data$i_man+data$i_don_ind+data$i_don_cor+data$i_sub
data$share <- (data$i_don_cor+data$i_don_ind)/data$total_income
data$share_individual <- (data$i_don_ind)/data$total_income
data$share_corporation <- (data$i_don_cor)/data$total_income
data$total_income_in_100k <- data$total_income/100000


# Full state names

data$state <- memisc::recode(as.factor(data$state),"Brandenburg" <- "BB","Berlin"<- "BE",
                             "Baden-Württemberg" <- "BW",
                             "Bavaria" <- "BY",
                             "Bremen" <- "HB",
                             "Hesse" <- "HE",
                             "Hamburg" <- "HH",
                             "Meckl.-Vorpommern" <- "MV",
                             "Lower Saxony" <- "NI",
                             "North Rhine-Westph." <- "NW",
                             "Rhineland-Palatinate" <- "RP",
                             "Schleswig-Holstein" <- "SH",
                             "Saarland" <- "SL",
                             "Saxony" <- "SN",
                             "Saxony-Anhalt" <- "ST",
                             "Thuringia" <- "TH")


data$Party <- data$party

# Build graphs

share_individual <- ggplot(data = data ,
                           aes(x = year,
                               y = share_individual,
                               shape=Party)) + theme_bw() + xlab("Year") + geom_point() + ylab("Share of Individual Donations in Total Income") + facet_wrap( ~ state, nrow = 3)  + scale_x_continuous(limits = c(2009, 2017), breaks = c(2009,2011,2013,2015,2017)) + geom_line()

# Plot

pdf(file = "Figure_A1.pdf",
    width = 10,
    height = 6,
    onefile=F)

share_individual + theme(axis.text.x = element_text(size = 6, angle = 45,hjust = 1.25, vjust = 1.5))

dev.off()

# Figure A2 ---------------------------------------------------------------

# Build graph

share_corporation <- ggplot(data = data ,
                            aes(x = year,
                                y = share_corporation,
                                shape=Party)) + theme_bw() + xlab("Year") + geom_point() + ylab("Share of Corporative Donations in Total Income") + facet_wrap( ~ state, nrow = 3)  + scale_x_continuous(limits = c(2009, 2017), breaks = c(2009,2011,2013,2015,2017)) + geom_line()

# Plot

pdf(file = "Figure_A2.pdf",
    width = 10,
    height = 6,
    onefile=F)

share_corporation + theme(axis.text.x = element_text(size = 6, angle = 45,hjust = 1.25, vjust = 1.5))

dev.off()

# Figure A3 ---------------------------------------------------------------

# Build graph

absolute_income <- ggplot(data = data ,
                          aes(x = year,
                              y = total_income_in_100k,
                              shape=Party)) + theme_bw() + xlab("Year") + geom_point() + ylab("Total Income in 100,000k Euro") + facet_wrap( ~ state, nrow = 3,scales="free_y")  + scale_x_continuous(limits = c(2009, 2017), breaks = c(2009,2011,2013,2015,2017)) + geom_line()

# Plot

pdf(file = "Figure_A3.pdf",
    width = 10,
    height = 6,
    onefile=F)

absolute_income + theme(axis.text.x = element_text(size = 6, angle = 45,hjust = 1.25, vjust = 1.5))

dev.off()

# Figure A4 ---------------------------------------------------------------

ENEP_data <- aggregate(ENP~year+state, data=data, FUN=mean)

# Build plot

ENEP <- ggplot(data = ENEP_data ,
                aes(x = year,
                    y = ENP,)) + theme_bw() + xlab("Year") + geom_point() + ylab("Effective Number of Electoral Parties") + facet_wrap( ~ state, nrow = 3)  + geom_hline(aes(yintercept =4.5)) + scale_x_continuous(limits = c(2009, 2017), breaks = c(2009,2011,2013,2015,2017)) 



pdf(file = "Figure_A4.pdf",
    width = 10,
    height = 6,
    onefile=F)

ENEP + theme(axis.text.x = element_text(size = 6, angle = 45,hjust = 1.25, vjust = 1.5))

dev.off()


# Figure A5 ---------------------------------------------------------------

# prepare data

data$i_don <-
  data$i_don_ind + data$i_don_cor # building overall donation variable
data$inh_1000 <-
  data$inh * 1000 # getting overall number of inhabitants
data$GDP_1000000 <-
  data$gdp * 1000000 # getting overall amount of GDP

# set variables to factors

data$fed_opp <-
  car::recode(data$federal_gov, "0=1;1=0", as.factor = TRUE)
data$state_gov <- as.factor(data$state_gov)
data$state_el <- as.factor(data$state_el)
data$federal_el <- as.factor(data$federal_el)

# round DVs (for negative binomial model)

data$i_don <- round(data$i_don, digits = 0)
data$i_don_ind <- round(data$i_don_ind, digits = 0)

data$i_don_cor <- round(data$i_don_cor, digits = 0)

# rescale IVs

data$inh_1000sd <-  scales::rescale(data$inh_1000,to = c(0,100))
data$inh_1000sd <- round(data$inh_1000sd, digits = 0)

data$GDP_1000000sd <-  scales::rescale(data$GDP_1000000,to = c(0,100))
data$GDP_1000000sd <- round(data$GDP_1000000sd, digits = 0)

# transform ENP

data$ENPx10 <-
  data$ENP * 10 # multiply ENP times 10 to receive integers
data$ENPx10 <- round(data$ENPx10, digits = 0)

# Run model

model_don_1_2 <-
  glmmTMB(
    i_don ~ ENPx10 + state_gov + fed_opp + state_el + federal_el + ENPx10 *
      state_el
    + GDP_1000000sd + inh_1000sd + (1 |
                                      state) + (1 |
                                                  party/state),
    family = nbinom2(link = "log"),
    REML = FALSE,
    data = data
  )

margeff_model_all <-
  ggeffects::ggpredict(model_don_1_2, c("ENPx10 [all]", "state_el"))
margeff_model_all_plot <-
  plot(margeff_model_all, colors = "bw", show.legend = TRUE) + labs(
    x = "ENEP*10",
    y = "All donations in euro",
    title = "",
    linetype  = "State election"
  ) + scale_linetype_manual(values = 15:17, labels = c("no","yes"))

# Plot

pdf(file = "Figure_A5.pdf",
    width = 5,
    height = 5,
    onefile=F)

margeff_model_all_plot

dev.off()

# Table A1 ----------------------------------------------------------------

data %>% dplyr::select(
  year,
  state_gov,
  fed_opp,
  state_el,
  gdp,
  inh,
  ENP,
  i_fees,
  i_don_ind,
  i_don_cor,
  i_man,
  i_sub
) %>% 
  mutate(across(where(is.factor), as.numeric)) -> 
  descriptives

stargazer(
  descriptives,
  type = "latex",
  summary = T,
  covariate.labels = c(
    "Year",
    "State Government",
    "Federal government",
    "State election",
    "GDP in Million euro",
    "Inhabitants in 1000 persons",
    "ENEP",
    "Income: Fees",
    "Income: Individual Donations",
    "Income: Corporative Donations",
    "Income: Mandate Holder Fees",
    "Income: State Subsidies"
  )
)

# Figure A6 ---------------------------------------------------------------

pdf(file = "Figure_A6.pdf",
    width = 8,
    height = 5,
    onefile=F)

ggplot(data, aes(x = i_don_ind)) + geom_histogram(bins = 500) + labs(x =
                                                                                  "Individual donations in euro", y = "Count") + theme_light()
dev.off()

# Figure A7 ---------------------------------------------------------------

pdf(file = "Figure_A7.pdf",
    width = 8,
    height = 5,
    onefile=F)

ggplot(data, aes(x = i_don_cor)) + geom_histogram(bins = 500) + labs(x =
                                                                                  "Corporate donations in euro", y = "Count") + theme_light()

dev.off()

